Report Only Show Code Show All

RRSG 2019 Challenge

Convention

BLUE BOX Steps for RRSG 2019 challenge tasks
  • Running cells followed by the blue boxes is enough to create RRSG 2019 challenge tasks.
YELLOW BOX Supplementary operations
  • Cells following yellow boxes contain supplementary implementations to generate interactive figures and to perform reconstruction using BART modules.
RED BOX Warnings to the user

Warnings to the user regarding the execution of the cells will appear in the red boxes:

  • The order in which the cells in a Jupyter Notebook are executed matters to ensure that all needed variables are in scope to run upcoming tasks.
  • Running some of the cells more than once may hamper the functionality of the remaining cells.

Clarifications

  1. How?

  2. This notebooks uses Script of Scripts (SoS) to exchange data between Octave and Python subkernels. This is achieved by the %get magic command placed at the very beginning of the Python3 cells. For example:

%get var_in_octave --from Octave

Following this %get magic command, var_in_octave variable will be also available to Python3 cells.

Initialize neccesary modules

INITIALIZE This notebook combines Octave and Python kernels for processing and visualization purposes, respectively. The following 4 cells are responsible for:
  • Loading Octave's image and optimization packages (included in the container, please see apt.txt)
  • Importing Python packages (included in the container, please see postBuild)
  • Defining two python functions to get a heatmap trace
  • Mex two .c implementations for calculating density compensation to use them in Octave
RUN FIRST The following 4 cells are run on opening to load neccesary modules!
In [3]:
%use Octave
disp('Executed automatically');
disp('Cell to load octave packages. To show, select the cell and click arrow up icon in the toolbar.');
pkg load image 
pkg load optim

Executed automatically
Cell to load octave packages. To show, select the cell and click arrow up icon in the toolbar.
In [4]:
# Load python packages
print('Executed automatically');
print('Cell to import python modules. To show, select the cell and click arrow up icon in the toolbar.');
import plotly.plotly as py
import plotly.graph_objs as go
import plotly_express as px
from plotly import tools
import numpy as np
import ipywidgets as widgets
import math
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=True)

Executed automatically
Cell to import python modules. To show, select the cell and click arrow up icon in the toolbar.
In [6]:
# This is a Python3 cell to create an interactive figure. 
# Here we use %get magic function to migrate variables from the Octave workspace into Python
print('Executed automatically');
print('Cell to define some Plotly functions. To show, select the cell and click arrow up icon in the toolbar.');

def heatmap_trace(inp, name, xlen, ylen, clrmax, clrmin, clrscale):
    trace = go.Heatmap( z =inp,
                        y = list(range(xlen-1)),
                        x = list(range(ylen-1)),
                        colorscale=clrscale,
                        showscale = True,
                        zmax=clrmax,
                        zmin=clrmin,
                        colorbar=dict(
                        tickfont=dict(
                            size=14,
                            color='white'
                        )),
                        name = name);
    return trace

def heatmap_trace2(inp, name, xlen, ylen, clrscale):
    trace = go.Heatmap( z =inp,
                        y = list(range(xlen-1)),
                        x = list(range(ylen-1)),
                        colorscale=clrscale,
                        showscale = False,
                        name = name);
    return trace

Executed automatically
Cell to define some Plotly functions. To show, select the cell and click arrow up icon in the toolbar.
In [7]:
% Mex c files for gridding by Brian Hargreaves and Philip Beatty 
% http://mrsrl.stanford.edu/~brian/gridding/
disp('Executed automatically');
disp('Cell to MEX c functions for dcf. To show, select the cell and click arrow up icon in the toolbar.');

mex gridlut_mex.c
mex calcdcflut_mex.c

Executed automatically
Cell to MEX c functions for dcf. To show, select the cell and click arrow up icon in the toolbar.

1. Load ISMRM RRSG 2019 challenge data.

1.1 Please see this OSF page for details.

Data is downloaded from the OSF and added to the portable software environment using postBuild configuration file.

In [11]:
load('/tmp/rrsg_challenge/brain_radial_96proj_12ch.mat');
whos % Show variables in the current scope

Variables in the current scope:

   Attr Name            Size                     Bytes  Class
   ==== ====            ====                     =====  ===== 
        ans             1x1                          8  double
   c    rawdata        12x96x512               4718592  single
        trajectory     96x512x3                 589824  single

Total is 737281 elements using 5308424 bytes
1.2 Change data order to follow BART's dimension convention
Warning Do not run the following cell more than once after loading data (previous cell).
Otherwise, data will be permutted once again and won't be following BART's convention anymore. --> You can refer to [this documentation](https://buildmedia.readthedocs.org/media/pdf/bart-doc/latest/bart-doc.pdf) for BART's dimension conventions.
In [12]:
rawdata = permute(rawdata,[4,3,2,1]); 
trajectory = permute(trajectory,[3,2,1]);
[~,nFE,nSpokes,nCh] = size(rawdata);
whos

Variables in the current scope:

   Attr Name            Size                     Bytes  Class
   ==== ====            ====                     =====  ===== 
        ans             1x1                          8  double
        nCh             1x1                          8  double
        nFE             1x1                          8  double
        nSpokes         1x1                          8  double
   c    rawdata         1x512x96x12            4718592  single
        trajectory      3x512x96                589824  single

Total is 737284 elements using 5308448 bytes
1.3 Display raw data from each channel using Pyhton.
In [13]:
%get rawdata --from Octave
%get nFE --from Octave 
%get nSpokes --from Octave
%get nCh --from Octave
print('Click on this cell and execute to pass variables from Octave to Python')
# This is a data exchange cell Octave --> Python3

Click on this cell and execute to pass variables from Octave to Python
In [15]:
# This cell plots rawdata from each channel
# The execution of this cell takes a bit longer than conventional 2D plots (e.g. matplotlib) as each datapoint is represented on a heatmap for interactivity using Plotly. 

rawdata = np.squeeze(rawdata)

fig = tools.make_subplots(rows=2, cols=6, print_grid=False, horizontal_spacing = 0.02, vertical_spacing = 0.02)

traces = []
iter = 0
for ii in range(2):
    for zz in range(6):
        cur_trace = heatmap_trace(np.log(1+abs(rawdata[:,:,iter])), 'Channel: '+ str(iter+1), nSpokes, nFE, 0.0001, -0, 'Viridis')
        fig.append_trace(cur_trace, ii+1, zz+1)
        iter += 1

fig['layout'].update(height=600, width=800, title=dict(text='<b>Raw data from each channel</b> <br> <i>Hover to see channel number and data points.</i>',font=dict(color="#ffc107")),paper_bgcolor='#000000')

for ii in range(12):
    exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels = False, ticks = \'\')')
    exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels = False, ticks = \'\')')

# plot(fig,filename='rawdata.html')
print('Execute this cell to regenerate the following figure.')
iplot(fig)

Execute this cell to regenerate the following figure.
In [53]:
from IPython.display import IFrame

IFrame(src='./rawdata.html', width=820, height=600,scrolling='no')
Out[53]:
1.4 Display k-space data on trajectories
In [31]:
# First we will parse data on Octave 

# Use bart to obtain root sum of squares of the rawdata over channels
rd = real(rawdata) + 1i*imag(rawdata);
rd = bart('rss 8', rd);
clr = squeeze(log(rd)); clear rd;
trajx = squeeze(trajectory(1,:,:));
trajy = squeeze(trajectory(2,:,:));
disp('Execute this cell before the follwing (if needed)');

Execute this cell before the follwing (if needed)
In [32]:
%get trajx --from Octave
%get trajy --from Octave
%get clr --from Octave

# The reason trajx, trajy and trajz are not vectorized in the Octave cell above is that, 
# currently SoS data exchange gets stuck while transferring large vectors. 

trajx = np.reshape(np.array(trajx),(nFE*nSpokes))
trajy = np.reshape(np.array(trajy),(nFE*nSpokes))
clr = np.reshape(np.array(clr),(nFE*nSpokes))

trac = go.Scatter(
    x  = trajx,
    y =  trajy,
    mode='markers',
    marker=dict(
        color = clr, #set color equal to a variable
        colorscale= 'Viridis',
        showscale=True,
        colorbar=dict(
        tickfont=dict(
            size=14,
            color='white'
        ))
    ),
    
)

layout = go.Layout(
    autosize=False,
    width=800,
    height=800,
    paper_bgcolor='#000000',
    plot_bgcolor='#000000',
    xaxis = go.layout.XAxis(
    gridcolor= '#283442',
         tickfont=dict(
            size=14,
            color='white'
        )
    ),
    yaxis = go.layout.YAxis(
    gridcolor= '#283442',
        tickfont=dict(
            size=14,
            color='white'
        )
    ),
    hovermode = 'closest',
    title=dict(text='<b>Raw data projected on k-space trajectory</b> <br> <i>Hover to see the locations of each sample</i>',font=dict(color="#ffc107"))
   
)

fig = go.Figure(data=[trac], layout=layout)

#plot(fig, filename='kspace.html')
print('Execute this cell (after previous one) to reproduce the following figure.')
iplot(fig)

Execute this cell (after previous one) to reproduce the following figure.
In [29]:
from IPython.display import IFrame

IFrame(src='./kspace.html', width=820, height=820,scrolling='no')
Out[29]:

2. Estimate coil sensitivities using BART

2.1 Use BART to estimate coil sensitivities
In [35]:
% Perform NUFFT to interpolate data onto cartesian grid. 
% -d denotes dimension (x:y:z, which is 300X300X1)
% -i sets the transform type to inverse
% -l enables L2 regularization  
% -t enables Toeplitz embedding
%  trajectory is the k-space locations of the acquired samples
%  rawdata is the sampled raw data

im = bart('nufft -d300:300:1 -i -l -t',trajectory,rawdata);

% Transform sensitivity maps to cartesian k-space using FFT
% -u denotes unitary transform
% 7 sets the bitmask level

im_ks = bart('fft -u 7', im);

% For details regarding BART's ECALIB, please see the implementation notes section.

calib = bart('ecalib -m1 -I',im_ks);

Done.
Done.
2.2 Display sensitivity profiles per channel
In [36]:
%get calib --from Octave
print('Transfer coil sensitivities to Python.')
# This is a data exchange cell. Octave --> Python3

Transfer coil sensitivities to Python.
In [37]:
calib = np.squeeze(calib)
fig = tools.make_subplots(rows=3, cols=4, print_grid=False, horizontal_spacing = 0.02, vertical_spacing = 0.02)

traces = []
iter = 0
for ii in range(3):
    for zz in range(4):
        mx = np.max(np.log(1+abs(calib[:,:,iter])))
        mn = np.min(np.log(1+abs(calib[:,:,iter])))
        cur_trace = heatmap_trace(np.log(1+abs(calib[:,:,iter])), 'Channel: '+ str(iter+1), 300, 300, mn,mx , 'Viridis')
        fig.append_trace(cur_trace, ii+1, zz+1)
        iter += 1

fig['layout'].update(height=750, width=800, title=dict(text='<b>Sensitivity profile of each channel</b> <br> <i>Hover to see channel number and data points.</i>',font=dict(color="#ffc107")),paper_bgcolor='#000000')

for ii in range(12):
    exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels=False,showline = False, ticks = \'\')')
    exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showline=False, showticklabels = False, ticks = \'\')')

# plot(fig,filename='calibrations.html')
print('Execute to renegerate coil sensitivity figure.');
iplot(fig)

Execute to renegerate coil sensitivity figure.
In [38]:
from IPython.display import IFrame

IFrame(src='./calibrations.html', width=820, height=820,scrolling='no')
Out[38]:

3. Reconstruction function

3.1 Use SENSE recon (both cartesian and non-cartesian) provided by BART for demonstration
3.1.1 Cartesian SENSE (BART) using k-space data obtained by fft of nufft
In [40]:
% Use gridded data for SENSE recon --> BART PICS 
bart_SENSE = bart('pics -l2', im_ks, calib);

Size: 90000 Samples: 90000 Acc: 1.00
l2 regularization: 0.000000
conjugate gradients
Total Time: 1.530570
3.1.2 Non-Cartesian SENSE (BART) using provided data
In [41]:
% Use non-cartesian SENSE recon --> BART PICS 
bart_SENSE2 = bart('pics -t',trajectory, rawdata, calib);

conjugate gradients
Total Time: 4.487198
3.1.3 Compare BART reconstructions from 3.0.1 and 3.0.2
RUN FIRST Depends on 3.1.1 and 3.1.2
In [42]:
%get bart_SENSE --from Octave 
%get bart_SENSE2 --from Octave

fig = tools.make_subplots(rows=1, cols=2, print_grid=False,vertical_spacing = 0.02)

fig.append_trace(heatmap_trace2(abs(bart_SENSE), 'BART pics -l2', 300, 300,'Greys'),1,1)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2), 'BART pics -t', 300, 300,'Greys'),1,2)

fig['layout'].update(height=400, width=600, title=dict(text='<b>Cartesian vs non-cartesian SENSE in BART</b>',font=dict(color="#ffc107")),paper_bgcolor='#000000')

for ii in range(2):
    exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels = False, ticks = \'\')')
    exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels = False, ticks = \'\')')

#plot(fig,filename='bart_cart_non.html')
print('Execute this cell to reproduce following figure.')
iplot(fig)

Execute this cell to reproduce following figure.
In [47]:
from IPython.display import IFrame

IFrame(src='./bart_cart_non.html', width=800, height=420,scrolling='no')
Out[47]:
3.2 Subsample provided data by factors of 2, 3 and 4.
In [54]:
% Please see HelperFunctions folder for subSample.m 

[outRD_x2, outTR_x2] = subSample(rawdata,trajectory,2,nSpokes);

[outRD_x3, outTR_x3] = subSample(rawdata,trajectory,3,nSpokes);

[outRD_x4, outTR_x4] = subSample(rawdata,trajectory,4,nSpokes);
3.3 Use non-cartesian SENSE recon (BART) to observe the effects of subsampling
3.3.1 Perform non-cartesian SENSE recon (BART)
In [55]:
% BART non-cartesian sense outputs 
bart_SENSE2_x2 = bart('pics -t',outTR_x2, outRD_x2, calib);
bart_SENSE2_x3 = bart('pics -t',outTR_x3, outRD_x3, calib);
bart_SENSE2_x4 = bart('pics -t',outTR_x4, outRD_x4, calib);

conjugate gradients
Total Time: 3.900349
conjugate gradients
Total Time: 4.888830
conjugate gradients
Total Time: 4.374562
3.3.2 Display: Compare reconstructed images
In [56]:
%get bart_SENSE --from Octave 
%get bart_SENSE2_x2 --from Octave
%get bart_SENSE2_x3 --from Octave
%get bart_SENSE2_x4 --from Octave
print('Execut to pass variables from Octave to Pyhton.')
# Data exchange cell. Octave --> Python3

Execut to pass variables from Octave to Pyhton.
In [57]:
fig = tools.make_subplots(rows=2, cols=2, print_grid=False,vertical_spacing = 0.02, horizontal_spacing=0.04)

fig.append_trace(heatmap_trace2(abs(bart_SENSE2), 'Original', 300, 300,'Greys'),1,1)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2_x2), '2X', 300, 300,'Greys'),1,2)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2_x3), '3X', 300, 300,'Greys'),2,1)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2_x4), '4X', 300, 300,'Greys'),2,2)

fig['layout'].update(height=600, width=600, title=dict(text='<b>Comparison of BART non-cartesian SENSE outputs <br> at different undersampling rates</b>',font=dict(color="#ffc107")),paper_bgcolor='#000000')

for ii in range(4):
    exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, showline=False, zeroline = False, showticklabels = False, ticks = \'\')')
    exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showline=False, showticklabels = False, ticks = \'\')')

#plot(fig,filename='sub_bart_compare.html')
print('Execute to reproduce the following figure');
iplot(fig)

Execute to reproduce the following figure
In [59]:
from IPython.display import IFrame

IFrame(src='./sub_bart_compare.html', width=620, height=620,scrolling='no')
Out[59]:
Section 4 Implementation of the reconstruction algorithm
  • To better relate implementation to the paper, reminders and implementation notes (following the original notation) are added.

  Reminder: Figure-1 from the paper 



Implementation of iterative image reconstruction.

Conjugate gradient (CG) iteration is controlled by the central CG process. It is initialized by MR data originating from N receiver channels (1,2,3,…N), acquired with an arbitrary k‐space trajectory.
Separately for each channel, these data undergo processing similar to conventional gridding reconstruction, i.e., sampling density correction (D), and resampling along a Cartesian grid, followed by FFT (FT1). The resulting images are individually multiplied by complex conjugate coil sensitivity and summed. After subsequent intensity correction (I), the sum image represents the vector a as defined by Eq. [25]. After initialization with a, the CG process iteratively calculates a progression of images, which converges towards exact reconstruction. For each iteration step, a current residuum image vector needs to be multiplied by the matrix IEH DEI. This is performed by the loop which starts from the CG box. After initial intensity correction (I), the processing is continued separately for each receiver coil. First, the intensity corrected residuum image is multiplied by individual coil sensitivity (Si). The results are transformed into k‐space by FFT and resampled along the experimental k‐space trajectory (FT2), resulting in a set of multiple‐coil k‐space data similar to that obtained experimentally. The following steps are equivalent to those carried out with the original data, yielding an intermediate image, which is fed back into the CG process. Here a refined approximate solution is calculated, which serves for further refinement by continued iteration.
As soon as the current approximation is sufficiently accurate, it is output and undergoes final intensity correction and k‐space filtering.

  Implementation notes  

This section relates the present implementation with the Figure-1 using the same notations.

  •  FT1 

    operation is implemented using BART's NUFFT:

This function transforms non-cartesian k-space data (rawdata and trajectory) into image domain (im):

im = bart('nufft -d x:y:z -i -l -t',trajectory,rawdata); 

% -d denotes dimensions in x:y:z
% -i sets the transform type to inverse
% -l enables L2 regularization  
% -t enables Toeplitz embedding
%  trajectory is the k-space locations of the acquired samples
%  rawdata is the sampled raw data

  •  S$\gamma$ 

    coil sensitivities are estimated using BART's ECALIB. See section 2.1

This function creates sensitivity maps (stored in the calib variable) from the raw k-space data on cartesian grid (sens_maps_ks):

calib = bart('ecalib -m1 -I',sens_maps_ks);

% -I enables intensity correction
% -m1 sets number of maps to one

  •  FT2 

    is implemented using BART's NUFFT:

This function transforms image im back to a non-uniform k-space of (radial, which is defined by trajectory variable in this case) nu_ks:

nu_ks = bart('nufft',traj,im);

  •  D 

    is the operation of scaling rawdata with density correction matrix.

To obtain density correction factor dcf pertaining to the rawdata, submodules (written in C) provided by Brian Hargreaves's gridding functions are used:

dcf = calcdcflut(trajectory,300); 

% trajectory is the k-space locations of the non-cartesian rawdata
% 300 is the recon matrix size (300x300)

  •  Sum 

    is implemented using BART's RRS:

This function simply performs sum of squares combination:

sum = bart('rss bitmask',multi_channel_input);

  •  F 

    is the k-space filtering, defined in 3.6
3.2 Calculate density compensation factor for operation

 D 

NOTE: This operation may take a few minutes.

A simpler method could have been applied here. Instead, I preferred using the submodules (written in C) provided by Brian Hargreaves's gridding functions to provide a good example of interoperability.

In [60]:
dcf = calcdcflut(trajectory,300); 
dcf = reshape(dcf,[3 nFE nSpokes]);

Warning: k-space trajectory should be complex-valued, or Nx2
3.3 Calculate intensity correction matrix for operation

 I 

In [61]:
# Root sum of square of sensitivities (estimated in section 2.1) from 12 channels 
I = abs(bart('rss 8',calib.*conj(calib)));
In [65]:
%get I --from Octave


trace  = heatmap_trace2(abs(I), 'I', 300, 300,'Viridis')
axis_template = dict(showgrid = False, zeroline = False,
             linecolor = 'black', showticklabels = False,
             ticks = '' )
layout = dict(height=500, width=500, title=dict(text='<b>Intensity correction matrix</b>',font=dict(color="skyblue")),paper_bgcolor='#000000',xaxis=axis_template,yaxis=axis_template)
figure = dict(data=[trace],layout=layout)


#plot(figure,filename='incor.html')
print('Execute to generate the following figure');

Execute to generate the following figure
In [68]:
from IPython.display import IFrame

IFrame(src='./incor.html', width=500, height=500,scrolling='no')
Out[68]:
3.4 Define operation

 E 

  • Function prototype
function E = opE(inp,S,traj,I)

  • Input arguments

    • inp: The image data (300,300)
    • S: Sensitivity profiles from 12 channels (300,300,12)
    • traj: k-space coordinates (3,512,96)
    • I: Intensity correction matrix

  • Output

    • E Please see the Figure-1 above

The following cell executes automatically to ensure that function definition exists.

In [70]:
function E = opE(inp,S,traj,I)

    inp = inp.*I; % Scale with the intensity correction matrix, as operation I precedes operation E (see Fig. 1)

    tmp = S.*inp; % Multiply the intensity corrected image with coil sensitivities. This will produce one image per channel stored in tmp variable

    E = bart('nufft',traj,tmp); % Transform back to the non-uniform k-space (see implementation notes for details)

end
3.5 Define operation

 EH 

  • Function prototype
function EH = opEH(dcf,inp,S,traj,I)

  • Input arguments

    • dcf: The density correction factor to scale rawdata
    • inp: Non-cartesian k-space samples from multiple channels (512,96,12)
    • S: Sensitivity profiles from 12 channels
    • traj: k-space coordinates (3,512,96)
    • I: Intensity correciton matrix

  • Output

    • EH Please see the Figure-1 above

The following cell executes automatically to ensure that function definition exists.

In [71]:
function EH = opEH(dcf,inp,S,traj,I)

  tmp = bart('nufft -d300:300:1 -i -l -t',traj,inp.*sqrt(dcf(1,:,:))); % NUFFT to image domain (see implementation notes)
  
  Sstar = conj(S(:,:,1,:)); %  Complex conjugate of sensitivity profiles (see Fig. 1)
  
  tmp2 = bart('rss 8 ',tmp.*Sstar); % Images scaled by complex conjugate of sensitivity profiles and SOS combined
  
  EH = tmp2.*I; % Scale output image by intensity correction

end
3.6 Define operation

 F 

The following cell executes automatically to ensure that function definition exists.

In [72]:
function out = opF(im,beta,r,N,I)
    
    im = im./I;
    %imagesc(im);
    im(isnan(im)==1) = 0;
    f_k = zeros(N,N);
    f_k(N/2+1,N/2+1) = 1;
    f_k = bwdist(f_k);
    f_k = 0.5 + 1/pi.*atan(beta.*((r-abs(f_k))/r));
    f_k = fftshift(f_k);
    im_k = fft2(im);

    out = abs(ifft2(ifftshift(im_k.*f_k)));

end
3.7 Define iterative algorithm

The following cell executes automatically to ensure that function definition exists.

In [73]:
function [b,deltas] = cg_solve(a,I,S,dcf,maxstep,trajectory)

    p = a;
    r = a;
    b = zeros(300,300);

    deltas = zeros(maxstep,1);
    for ii=1:maxstep

    disp(['Iteration -->' num2str(ii)]);

    delta = r(:)'*r(:)/(a(:)'*a(:));
    
    disp(delta);
    deltas(ii) = delta;

    E = opE(p,S,trajectory);

    q = opEH(dcf,E,S,trajectory,I);

    % dot(r,r) is equivalent to r(:)'*r(:). Used dot for easy reading. 

    term = dot(r,r)/dot(p,q);

    b = b + term*p;

    rprev = r;

    r = r - term*q;

    term2 = dot(r,r)/dot(rprev,rprev);

    p = r + term2*p;

    end

    
end
3.8 Define main call to the iterative solution

The following cell executes automatically to ensure that function definition exists.

In [74]:
function [im,deltas] = main_sense(rawdata, calib, trajectory,I,dcf,maxiter)

    a = opEH(dcf,rawdata,calib,trajectory,I);
    [im,deltas] = cg_solve(a,I,calib,dcf,maxiter,trajectory);
    im = opF(im,100,40,300,I);

end
4.1 Perform recon with different subsamplings
In [75]:
[im,deltas] = main_sense(rawdata, calib, trajectory,I,dcf,10);
[im_x2,deltas_x2] = main_sense(outRD_x2, calib, outTR_x2,I,dcf(:,:,1:2:end),10);
[im_x3,deltas_x3] = main_sense(outRD_x3, calib, outTR_x3,I,dcf(:,:,1:3:end),10);
[im_x4,deltas_x4] = main_sense(outRD_x4, calib, outTR_x4,I,dcf(:,:,1:4:end),10);

Done.
Iteration -->1
 1
Done.
Done.
Iteration -->2
 0.17712
Done.
Done.
Iteration -->3
 0.13987
Done.
Done.
Iteration -->4
 0.098370
Done.
Done.
Iteration -->5
 0.096087
Done.
Done.
Iteration -->6
 0.096355
Done.
Done.
Iteration -->7
 0.096350
Done.
Done.
Iteration -->8
 0.096327
Done.
Done.
Iteration -->9
 0.096307
Done.
Done.
Iteration -->10
 0.096294
Done.
Done.
Done.
Iteration -->1
 1
Done.
Done.
Iteration -->2
 0.17666
Done.
Done.
Iteration -->3
 0.13934
Done.
Done.
Iteration -->4
 0.098537
Done.
Done.
Iteration -->5
 0.096187
Done.
Done.
Iteration -->6
 0.096440
Done.
Done.
Iteration -->7
 0.096423
Done.
Done.
Iteration -->8
 0.096398
Done.
Done.
Iteration -->9
 0.096377
Done.
Done.
Iteration -->10
 0.096361
Done.
Done.
Done.
Iteration -->1
 1
Done.
Done.
Iteration -->2
 0.17628
Done.
Done.
Iteration -->3
 0.14899
Done.
Done.
Iteration -->4
 0.10486
Done.
Done.
Iteration -->5
 0.10359
Done.
Done.
Iteration -->6
 0.10341
Done.
Done.
Iteration -->7
 0.10333
Done.
Done.
Iteration -->8
 0.10329
Done.
Done.
Iteration -->9
 0.10327
Done.
Done.
Iteration -->10
 0.10326
Done.
Done.
Done.
Iteration -->1
 1
Done.
Done.
Iteration -->2
 0.17516
Done.
Done.
Iteration -->3
 0.15006
Done.
Done.
Iteration -->4
 0.10419
Done.
Done.
Iteration -->5
 0.10327
Done.
Done.
Iteration -->6
 0.10296
Done.
Done.
Iteration -->7
 0.10284
Done.
Done.
Iteration -->8
 0.10279
Done.
Done.
Iteration -->9
 0.10277
Done.
Done.
Iteration -->10
 0.10276
Done.
Done.
In [82]:
%get deltas --from Octave
%get deltas_x2 --from Octave
%get deltas_x3 --from Octave
%get deltas_x4 --from Octave
In [89]:
iterations = [1,2,3,4,5,6,7,8,9,10];

sub2_1 = go.Scatter(
    x = iterations,
    y = np.log10(deltas),
    name = 'R=1',
    text = 'R=1',
    hoverinfo = 'x+y+text',
    line=dict(color='rgb(0,114,189)'),
)

sub2_2 = go.Scatter(
    x = iterations,
    y = np.log10(deltas_x2),
    name = 'R=2',
    text = 'R=2',
    hoverinfo = 'x+y+text',
    line=dict(color='rgb(217,83,25)'),
)

sub2_3 = go.Scatter(
    x = iterations,
    y = np.log10(deltas_x3),
    name = 'R=3',
    text = 'R=3',
    hoverinfo = 'x+y+text',
    line=dict(color='rgb(237,177,32)'),
)

sub2_4 = go.Scatter(
    x = iterations,
    y = np.log10(deltas_x4),
    name = 'R=4',
    text = 'R=4',
    hoverinfo = 'x+y+text',
    line=dict(color='rgb(126,47,142)'),
)

fig = tools.make_subplots(rows=1, cols=2, print_grid=False)
fig.append_trace(sub2_1, 1, 1)
fig.append_trace(sub2_2, 1, 1)
fig.append_trace(sub2_3, 1, 1)
fig.append_trace(sub2_4, 1, 1)



iplot(fig)
In [91]:
%get im --from Octave
%get im_x2 --from Octave
%get im_x3 --from Octave
%get im_x4 --from Octave


fig = tools.make_subplots(rows=2, cols=2, print_grid=False,vertical_spacing = 0.02, horizontal_spacing=0.04)

fig.append_trace(heatmap_trace2(abs(im), 'Original', 300, 300,'Greys'),1,1)
fig.append_trace(heatmap_trace2(abs(im_x2), '2X', 300, 300,'Greys'),1,2)
fig.append_trace(heatmap_trace2(abs(im_x3), '3X', 300, 300,'Greys'),2,1)
fig.append_trace(heatmap_trace2(abs(im_x4), '4X', 300, 300,'Greys'),2,2)

fig['layout'].update(height=600, width=600, title=dict(text='<b>Comparison of iterative algorithm <br> at different undersampling rates</b>',font=dict(color="#ffc107")),paper_bgcolor='#000000')

for ii in range(4):
    exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, showline=False, zeroline = False, showticklabels = False, ticks = \'\')')
    exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showline=False, showticklabels = False, ticks = \'\')')

#plot(fig,filename='sub_bart_compare.html')
iplot(fig)